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Abstract. We perform 2D numerical simulations of a 
magnetorotational explosion of a rotating magnetized gas 
cloud. We found that amplification of a toroidal magnetic 
field due to the differential rotation leads to a transforma- 
tion of the part of the rotational energy of the cloud to 
the radial kinetic energy. Simulations have been made for 
3 initial values of £ (the relation of magnetic energy to the 
gravitational energy of the cloud): £ = 1CP 2 , 1CP 4 , 1CP 6 . 
Part of the matter - ~ 7% of the mass of the cloud (~ 3.3% 
of the final gravitational energy of the cloud) - gets radial 
kinetic energy which is larger than its potential energy 
and can be thrown away to the infinity. It carries about 
30% of the initial angular momentum of the cloud. This 
effect is important for angular momentum loss in the pro- 
cesses of stellar formation, and for the magnetorotational 
mechanism of explosion suggested for supernovae. Simu- 
lations have been made on the basis of the Lagrangian 2D 
numerical implicit scheme on a triangular grid with grid 
reconstruction. 
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1. Introduction 



stars: supernovae ■ 



Magnetic field plays an important role in life of a star, 
especially at its birth and its death. A birth of a star in a 
galactic disk during the collapse of an interstellar cloud is 
complicated by large angular momentum of matter which 
will prevent formation of the star if it is not lost. The 
most realistic mechanism of the angular momentum loss 
is connected with a magnetic field which, twisting dur- 
ing differential rotation, leads to a flux of angular mo- 
mentum outside, pcrmiting a collapse of the main body 
(Bisnovatyi-Kogan et al. 1973, Ardeljan et al. 1996b). 

In the late stages of evolution of massive stars, loss of 
hydrodynamical stability initiates a collapse, which is fin- 



ished by supernovae explosions after formation of a stable 
neutron star. Investigation of supernovae explosions in dif- 
ferent nonmagnetized models revealed serious problems in 
transformation of the gravitational energy into the energy 
of explosion. In this situation magnetorotational mecha- 
nism of supernovae explosion, suggested by Bisnovatyi- 
Kogan (1970) could be an explanation of this phenomena. 

Magnetorotational phenomena in stellar envelopes are 
considered to be a main mechanism of an angular momen- 
tum loss from stars, and magnetized solar wind is a reason 
for a very slow solar rotation. 

2D calculations of collapse of a rotating magnetized 
star with a somewhat unrealistic magnetic field configu- 
ration have been performed by Le Blank & Wilson (1970) 
and Ohnishi (1983). The geometry of the outburst which 
they obtained, was different, but the energy of the out- 
burst was substantial in both cases. 

We investigate here magnetorotational phenomena in 
the collapsing rotating magnetized gas cloud initially with 
uniform density and angular momentum, which may be 
considered as a model of star formation. After mag- 
netorotational explosion a quasistationary, almost uni- 
formly slow rotating configuration is formed. The loss 
of angular momentum is connected with an outburst 
of the matter, carrying an energy which is ~ 3.3% of 
the final gravitational energy of the cloud. Our simula- 
tions show that amounts of ejected mass and energy are 
weakly dependent on the initial value of magnetic en- 
ergy (parameter £). The main difference in the results 
of 2D simulations for the different initial values of £ is 
in the duration of the process. This conclusion confirms 



results of ID simulations ( Bisnovatyi-Kogan et al. 1976 
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Ardeljan et al., 197£), with the time of the process ap- 
proximately proportional to the £ -1 ^ 2 . Similar processes 
related to a different initial model and equation of state 
are expected to act during supernovae explosions. 

Magnetorotational explosion of a rotating magnetized 
gas cloud in a 2D approac h has been investigated earlier 
by Ardeljan et al. (1996b)| for a divergency- free, but not 
force-free, magnetic field configuration. Magnetic forces 
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have not been balanced by a pressure distribution and 
gravitational forces. The effect of ejection of the part of 
the mass of the cloud has been found there, but use of an 
unbalanced initial field leads to artificial effects in a MHD 
flow and does not allow us to extend calculations as far as 
is needed. 

In this paper an initial configuration with magnetic 
forces balanced by pressure and gravitational forces dis- 
tribution has been constructed, allowing us to make calcu- 
lations for different sets of initial parameters and to extend 
them to the larger times. 

2. Basic equations 

Consider a set of magnetohydrodynamical equations with 
selfgravitation and with infinite conductivity (Landau & 
Lifshitz, 1984, |Bisnovatyi-Kogan, 1989| ): 



dx dp 

_ =U , -+pd 1VU : 
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dt 
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where ^ = Jj + u • V is the total time derivative, 
x = (r, tp, z), u is velocity vector, p is density, p is pres- 
sure, H = (H r , H v , H z ) is magnetic field vector, $ is 
gravitational potential, e is internal energy, G is gravita- 
tional constant, -ft is universal gas constant, 7 is adiabatic 
index, H <£> H is tensor of rank 2. 

Axial symmetry ( ) and symmetry to the equatorial 
plane [z = 0) are assumed. 

The problem is solved in the restricted domain. Out- 
side the domain, the density of the matter is zero, but 
poloidal components of magnetic field H r , H z can be non- 
zero. 

To write the set of equations in dimensionless form we 
choose the following scale values: 



Po 



= 1.492 • l(T 17 g/cm 3 , r = z = x = 3.81 • 10 16 cm, 



r = rxo, z — zxo, u — uuq, uq = yj inGpox^, 

p = ppo, e = ee , T = TT , $ = $$ , $ = ^Gp xl, 

, x 2 2.-2 rp 

to = — , Po = Pqu = poX Q t , lo - 

2 2.-2 tt 1 1 — 1 1/2 

e = u = x t , H = y/p^ = x t p Q ' . 
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Here the values with index zero are the scale factors and 
the functions under a tilde are dimensionless functions. 
The set of equations ([!]) can be written in the following 
nondimcntional form (the tilde being omitted): 
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Taking into account symmetry assumptions (^ = 0), 



the divergency of the tensor H 1 
the following form: 



H can be presented in 
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3. Collapse of a rotating magnetized gas cloud 

3.1. Formulation of the problem 

Consider a magnetized rotating gas cloud which is de- 
scribed by the set of equations (|l]) . All graphs and figures 
below are in a nondimcnsional form. At the initial moment 
(t = 0) the cloud is a rigidly rotating uniform gas sphere 

•"' - 1 " - 17 g/cm 3 , 7] = 5/3, 



(Fig. |l|) with the following parameters: 
r = 3.81 • 10 16 cm, p = 1.492 • 10" 



M = 1.73M = 3.457 • 10 33 g, u r = u z = 0, 



f3 r0 = ^ = , )4 . 
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Here the subscript "0" corresponds to the initial moment 
(t = 0), subscript "1" to the moment of the beginning of 
the evolution of magnetic field (t = t\ > 0). 

The assumptions about the symmetry of the problem 
are the same as in Section |^. 

The process can be divided into the following three 
qualitatively different stages. First is a pretty short hy- 
drodynamical collapse stage. At this stage the influence 
of the magnetic field on the process of the collapse of the 




Fig. 1. Grid at initial time t = 0. 



cloud can be neglected, because the initial poloidal mag- 
netic field is weak. In the short time of collapse the toroidal 
component of the magnetic field, which appears to be due 
to the arising differential rotation is also weak at this ini- 
tial stage. The second stage is the stage of a rather "long" 
twisting of magnetic field due to the differential rotation of 
the cloud. The final third stage starts with the appearance 
of a compression wave, moving from the inner parts of the 
cloud to its periphery along a steeply decreasing density 
background. Soon after its appearance, it transforms into 
the MHD shock wave, which can push out a light envelope 
of the protostar. Similar ejection following formation of a 
rapidly rotating neutron star and a differentially rotating 
envelope, can be interpreted as supernova explosion. 

Simulations for different £ consist of the calculation of 
oscillations of the cloud until formation of the differen- 
tially rotating equilibrium (without magnetic field), and 
subsequent inclusion and twisting of the magnetic field. 

3.2. Initial magnetic field configuration 

The initial magnetic field must satisfy the condition of 
absence of magnetic charges divH = 0. It also has to cor- 
respond to the boundary conditions of the problem. 

The best choice would be dipole or quadrupole. While 
they satisfy initial and boundary conditions, they have 
singularities in the origin of coordinates (r = 0, z = 0). 
Using such magnetic fields in numerical simulations can 
lead to loss of accuracy of calculations. 



To define the initial magnetic field we use the following 
method. We defined toroidal current j v in the central part 
of the core of the collapsed cloud by a formula: 
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After getting the defferentially rotating stationary so- 
lution for nonmagnetized cloud we use Bio-Savara law 
( A.l),( A.2),( A. 5) for calculation of the poloidal compo- 
nents of the magnetic field H r o, H z q (Fig. |^). This mag- 
netic field is divergency free, but it is not force-free and 
should be balanced at the initial moment. Then we use the 
following method: we "turn on" the poloidal magnetic field 
H r o, H z q, but "switch off" the equation for the evolution 
of the toroidal component H v in (jl|). Actually it means 
that we define H v = 0, —rr- = 0. From the physical point 
of view it means that we allow magnetic field lines to slip 
through the matter of the cloud in the toroidal direction. 
After "turning on " such a field, we let the cloud come 
to the steady state, where magnetic forces connected with 
the purely poloidal field are balanced by other forces. 

The calculated balanced configuration has the mag- 
netic field of quadrupole- like symmetry. For testing we run 
our code with this purely poloidal field for a large number 
of time steps (~ 10 3 ), during which the parameters of the 
cloud did not change. 
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Fig. 4. Density field at t x = 18.45. 



4. Results 

We describe results of simulations of the problem for initial 
£(ti) = 10~ 2 . Results of simulations for £,{t x ) = 10~ 4 , 1CT 6 
are qualitatively similar to the first case. The amounts 
of the ejected mass and energy are the same, ~ 7% of 
mass and ~ 3.3% of energy. The main difference betveen 
these three variants is the duration of the process. The 
lower the initial magnetic energy, the longer the evolution 
time up to the magnetorotational explosion. This feature 
of the results of the simulations in a 2D case is similar 
to the results of ID simulations of the magnetorotational 
mechanism for supernovae (Ardeljan, et.al. 1979). 
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Angular velocity distribution along r axis at ti 



We start calculations (t = 0) with the set of initial 
parameters (||)-(|^) (without magnetic field at the begin- 
ning). Our triangular grid at t = consists of 2200 knots 
and 4400 cells. During the calculations, the number of 
knots and cells deviates no more than 5% from its ini- 
tial values due to specially developed procedure (Ardel- 
jan, et.al. 1996a). After a number of oscillations the cloud 
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at t\ = 18.45 comes to the steady differentially rotating 
configuration. Lagrangian grid, density field and angular 
velocity levels are shown in Figs. [| ||, || Distribution of 
the density and angular velocity for the t\ = 18.45 along 
r-axis are given in Figs. ||, fj]. 

At t\ the cloud consists of a rapidly rotating dense core 
and a slowly rotating prolate envelope (Figs. |[ f7|). The 
radial part of kinetic energy of the cloud is less than 1.5% 
of its gravitational energy. The cloud can stay an infinite 
time in such differentially rotating stationary condition 
(if we neglect dissipation processes) . However inclusion of 
even a weak magnetic field leads to an essential change in 
its state. 

At t\ we " turn on" a poloidal magnetic field calculated 
as described above. Simultaneously at t\ we "switch off" 
the equation for the evolution of the toroidal magnetic 
field for a short time (as discussed in the section "Ini- 
tial magnetic field configuration" ) . After a few weak oscil- 
lations around equilibrium, the cloud comes to the steady 
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evolution of the toroidal magnetic field. The black colour 
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state configuration with a balanced magnetic force. These 
oscillations are clearly seen in the plot for the time evolu- 
tion of the poloidal magnetic energy (Fig. ||) This purely 
poloidal field becomes balanced at — 26.86 . At this 
time we "switch on" the equation for the toroidal com- 
ponent of magnetic field and let the magnetic force lines 
twist due to differential rotation. 

Due to the differential rotation of the cloud and very 
high (infinite) conductivity soon after t — t2, the toroidal 
component of magnetic field appears and grows with time. 
Taking into account the quadrupole-like symmetry of the 
initial poloidal magnetic fields the toroidal component ap- 
pears and has 2 extremes (Fig. ||). A maximum at the 
equatorial plane and a minimum at the periphery of the 
core of the cloud, close to the z-axis. These extrema at the 
initial stage of the evolution of the toroidal magnetic field 
correspond approximately to the extremes of the following 
scalar product: H • grad-^f, because the cloud at the mo- 
ment of " switching on" the equation for toroidal magnetic 
field was in a steady state condition and only this term in 
the equation for the H v determines the evolution of the 
toroidal field. 

Toroidal magnetic energy grows almost linearly start- 
ing from t2 up to reaching its maximum value (Fig. ^) after 
^2.5 rotations of the central core, which rotates almost 
rigidly. Increase of the toroidal magnetic field leads to sub- 
traction of angular momentum from the central parts of 
the cloud and additional contraction of the core of the 
cloud, meanwhile the envelope of the cloud starts to blow 
up. Contraction wave moving outwards appears at the pe- 
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Fig. 12. Toroidal magnetic field H v distribution along the 
r axis at t x = 28.886. 



riphery of the core of the cloud. This wave propagates 
along a decreasing density profile, increasing its amplitude 
and transforms to the MHD shock wave. 

At t — 28.0 a small outer part of the envelope of the 
cloud reaches the poloidal (radial) velocity (u r , v z ), which 
corresponds to a kinetic energy larger then its potential 
energy and can fly away to infinity. At this time maximal 
density in the center of the cloud is p m ax = 113.0 (Fig. |lo|) . 
The amount of matter thrown away by the MHD shock, 
grows with time. In figs. p^lJlq one can see the evolution 
of the velocity field and development of the process of 
ejection of the matter of the envelope of the cloud with 
time. 

The growth of the toroidal magnetic field leads to the 
formation of MHD shocks. To analyze the structure of 
the evolution of these shocks we consider a flow picture 
at the equatorial plane (along r-axis). The toroidal field 
at the r-axis, has the maximum at the periphery of the 
core of the cloud, approximately coinciding with the point 
where the radial velocity v r changes its sign (Fig. |TT| ) . 
At r « 0.15 (Fig. O) MHD shock is formed, and it is a 
slow MHD shock, because its velocity is lower than the 
local Alfvenic velocity. Angular velocity v v /r and H v is 
decreasing, when passing through this shock. The density 
p and temperature T are increasing. 



The matter which is on the right part of the maxi- 
mum of H v is moving outwards. The MHD shock, which 
is formed on the right side of the maximum of the 
toroidal field at the equatorial plane is a fast MHD 
shock wave, and its velocity is bigger than fast mag- 
netic sound speed before the front of the shock. The 
velocity of the MHD shock wave is smaller than the 
Alfvenic sound speed in the gas behind the shock, and 
bigger than the slow magnetic sound speed after the 



shock Kulikovskii & Ljubimov, 1962. The toroidal mag- 
netic field (Fig. |12[) and angular velocity (Fig. |l3|) grow 
behind the shock. Due to the transformation of the rota- 
tional energy into the kinetic energy of the radial motion, 
the angular velocity of the core of the cloud decreases with 
time (Fig.pl). At t = 28.886 the shock is at r w 1.25 (see 
for example Figjl^) . 

Similar slow MHD shock has been found in ID sim- 
ulations of magnetrotational supernova explosion made 
by Ardeljan, et.al. (1979), while the formation of the fast 
MHD shock was not presented. 

The time evolution of the parameters 



Mt) . m - *.<«> 



^"mag (^) 



(8) 



|£ gr (*)|— \E p (t)\'"-' \E p (t)\' 

for initial = 10 -2 , is given in Fig. (|l6|) . 

In Figs. [l7] and p[ time evolution is presented of the 
ejected mass of the envelope of the cloud, and the energy 
it carries away. Just after the beginning of the ejection the 
amount of the ejected mass grows approximately linearly 
with time. At the moment t — 33.57, its growth stops and 
until the end of the calculations these values (ejected mass 
and energy) do not change significantly. The amount of the 
ejected matter after this time is about 7% and it carries 
approximately 3.3% of the total energy of the cloud. At 
the moment of the evolution of the magnetorotational ex- 
plosion the rotation of the cloud changes essentially. The 
core of the cloud now rotates with the small angular ve- 
locity uj core ~ 0.85. The ejecting matter now looks like an 
expanding shell. At t = 28.0 (the moment of the "switch- 
ing on" of the evolution of the toroidal magnetic field, the 
"inner" (core and part of the envelope) 50% of the mass 
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of the cloud contained 26.8% of the angular momentum, 
at t = 33.57 the "inner" 50% of the mass of the cloud con- 
tained 7.2% of the angular momentum. At t = 28.0 the 
50% of the angular momentum are contained in the outer 
part of the envelope of the cloud (31.3% of the mass of the 
cloud). At t = 33.57 the 50% of the angular momentum 
is concentrated in the outer part of the envelope of the 
cloud (9.2% of the mass of the cloud). At the last plot of 
Fig. |l5| the structure of the cloud at the advanced stage 
of the magnetorotational explosion is given. At this stage 
the cloud consists of an outflowing envelope, a transitional 
region and an almost rigidly rotating core. 



In the paper by Ardeljan et al., 1979 it was found 
in ID calculations of the magnetorotational supernova 
explosion that during the evolution of the toroidal 
component and angular momentum transfer outwards, 
the central core of the cloud starts to rotate in 
the opposite direction to the initial one. Such mag- 
netorotational oscillations have been investigated ana- 



lytically by Bisnovatyi-Kogan et al. 1976 With respect 
to the star formation problems in slightly differ- 
ent formulation, these oscillations have been investi- 



gated in the papers by Moushovias fc Paleologou, 1979 



Moushovias fc Paleologou, 1980| . In our simulations we do 
not get the opposite rotation of the core of the cloud in 
the opposite direction. However the core loses significant 
part of its angular momentum due to magnetic breaking. 

5. Discussion 

The first 2D calculations of the magnetorotational explo- 
sion were performed by Le Blank & Wilson (1970), where 
matter was expelled in the form of jets. Such geometry 
of the outburst was the result of a specific rather unre- 
alistic choice of the magnetic field configuration, which 
was produced by a current ring, at an equator, out of a 
stellar centre, where the matter density was an order of 
magnitude less than a central one. The magnetic field of 
this ring had a zero radial component in the equatorial 
plane, and a magnetic pressure gradient was formed in 
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Fig. 18. Time evolution (in % to the total initial energy 
of the cloud) of the energy of the ejected matter of the 
envelope of the cloud. 



the z-direction due to the differential rotation, which al- 
most repeated the magnetic pressure gradient of the initial 
configuration. Such magnetic field structure had led to the 
matter stream pattern appearing preferentially along the 
symmetry axis of the magnetic field. 

In our calculations the magnetic field was created by 
the oppositely directed current rings in the central region 
of the star. The magnetic field had a quadrupole-like con- 
figuration with a maximum value close to the centre of the 
star, and large radial component at the equator. 

In the papers devoted to the simulations of the star 
formation problems (e.g. Ciolek & Mouschovias, 1995; 
Tomisaka, 1998) it is often supposed that the initial 
poloidal magnetic field has only an H z component. It is 
often suggested for the collapse problem of protostellar 
clouds. We should point out that the initial magnetic field 
of the configuration H$ z ^ 0, Ho r = or dipole-like mag- 
netic fields (Ho r = at the equatorial plane) can lead 
after, the evolution of the toroidal component, to the jet- 
like outflows directed presumably along the z-axis. How- 
ever application of the quadrupole-like initial poloidal field 
(H z = at the equatorial plane) can give radial ejections, 
as is seen in our results. 
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The initial magnetic field configuration in the rotat- 
ing gas cloud as a progenitor of the star formation may 
be rather complicated because of noncoincidence of mag- 
netic and rotational axes and possible action of the dy- 
namo mechanism. Due to this fact the prevailing of a 
quadrupole-like component in the initial model it seems 
to be possible in some cases also in the cloud. 

Observation of the solar corona during sunspot mini- 
mum has shown a north-south asymmetry, interpreted by 
Osherovich et al. (1984) as the presence of a significant 
magnetic quadrupole. Due to more rapid central increase 
of the quadrupole component (~ l/?" 4 ) in comparison with 
the dipole (~ 1/r 3 ), we may expect that quadrupole com- 
ponent may be larger than the dipole one in the central 
part of the star. Besides, we believe that the maximum 
of the magnetic field strength can not be situated at such 
large distance from the centre, as it was taken by Le Blank 
& Wilson (1970). Another 2D calculation, performed by 
Ohnishi (1983) for a radial initial magnetic field gave an 
outburst mainly in the equatorial plane, similar to our 
results. 

The amount of energy of the matter ejected during the 
magnetorotational explosion is sufficient for producing a 
supernova explosion after formation of the rapidly rotating 
neutron star with a differentially rotating envelope. Such 
calculations are now in progress. 

Appendix A: Numerical method 

We have used a numerical technique, based on a 
conservative (in the absence of gravitation), implicit 
first order of accuracy in space and time oper- 
ator difference scheme on a triangular grid with 
a grid reconstruction, developed and described in 



Ardeljan et al. (1987), Ardeljan fc Kosmachevskii (1995) 



This numerical method has been used in the paper by 
Ardeljan et al. (1996a)| for the investigation of the collapse 
of the nonmagnetized rotating cloud. 

For a numerical solution of the problem of the mag- 
netorotational explosion, we introduce a triangular grid, 
covering the restricted domain in r, z coordinates. 

We suppose that components of the velocity vector u 
and gravitational potential $ will be defined in all knots 
of the grid. The density p, pressure p, components of mag- 
netic field vector H will be defined in the cells and in the 
boundary knots of the grid. 

For the numerical simulation of the system of gravi- 
tational MHD equations the following method has been 
used. Instead of differential operators (div, grad, rot) we 
introduce their finite difference analogues. On the base of 
such operators a completely conservative scheme has been 
constructed. The scheme is implicit for all velocity com- 
ponents v r ,Vtp,v z and for the toroidal component of the 
magnetic field H v . The scheme is explicit for the poloidal 
magnetic field H r , H z . The explicitness of the scheme for 
H r ,H z does not introduce strong restriction on the time 



step, because during the evolution of the magnetic field its 
poloidal values do not change strongly, while the toroidal 
component appears and increases significantly with time. 
The scheme is explicit for the gravitational potential, but 
it was shown in Ardeljan et al. (1987)| that this explicit- 
ness does not introduce significant restrictions on the time 
step. 

A.l. Calculation of the initial and boundary values of the 
magnetic field 

Initial values of the poloidal components of the magnetic 
field H TQ , H Zo in the computational domain and its bound- 
ary values at the outer boundary H rq , H zq are calculated 
using the Bio-Savara law: 



H(M ) = i / 3 -^^dV M , 

c J Km AI a 
V 



(A.l) 



where J = j r e r +j ip e ip +j z e z - current density, e r , e^, e z - 
unit vectors of cylindrical system of coordinates, c - light 
velocity, Rmm - radius vector connecting Mq and M 
points, Rmm - length of the vector Rmm - 
In the cylindrical coordinates we have: 



H(M ) 



1 



rds 



J x R 



Rmm 



°-dp 



(A.2) 



where S is the computational domain in r, z coordinates. 
On the triangular grid the integral J f(r, z)rds can be 



approximated by the following sum: 
/ f{r,z)rds& ^ nsifir^Zi) 



(A.3) 



where u>i is a set of triangular cells, is r coordinate of 
the center of the cell i, Zi is z coordinate of the center of 
the cell i, s, - area of the grid cell i. 

For the integration of the vector- function in cylindrical 
coordinates we use the following formula: 



f2 



f r (p)e r + f v {ip)e v + f z (p)e z 



dip 



= e r (p 2 ) J fr(f) cos(ip 2 - <p) + ftp(v) sin(p 2 - p) dip 

vi 



¥2 

+ e v ((f2) J f v (p) cos(p 2 -<p)- f r (p) sin(p 2 ~ p) 

Vl 

V2 



dip 



+ e z (p 2 ) J fz(p)dp. 



(A.4) 
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Using (A.4) for (A.2), after transformations we get: 



A. 2. Testing 



H r (r, z) = ^2 r i s i j ip i< 



(z - Zi)- 
I 

+ (z + Zi)^z 

H v (r,z)=0 
H z (r,z) 



2-k 



^U(-k 2 ,h, 



F(k u 



^ ft 2 



^2 r i s d<pi\n- s 'n.{- k 2 ,h,^) 



+ r- 



k 2 n ( ~~ 



fi-gllf - fcf, fc 2 , ~) 



7T, 

2' 



Testing of the numerical technique, applied in this paper, 
without magnetic field has been described in detail by 



Ardeljan et al. (1996a 



In MHD case we add the following tests: 

1. We calculated a spherically symmetrical collapse of a 
nonrotating gas cloud without a magnetic field until 
the cloud comes to the steady state condition. Af- 
ter that we "turned on" force- free and divergency- free 
, magnetic field H r = 0, H v = 0, H z = const. This mag- 

netic field does not change the condition of the cloud 
after hundreds of numerical time steps. 
(A. 5) 2. The piston problem^ Consider the flat piston, which 
is pushed along the r-axis (our domain for this test 
was a rectangle with the following coordinates: (ro = 
1000, zq = 0; n = 1001, z\ = 1), in this case we 
can consider r, z as Cartesian coordinates with high 
accuracy) into the gas with following characteristics: 
infinite conductivity [a = oo), equation of state: 

p = 5RpT, e = 9fJT/( 7 - 1), 5K = 1, 7 = 5/3. 
At t = the gas has the following parameters: 



2-k 



^(-klk 2 ,-] 



Here 



, ai 2 = {n + r) 2 + (z- Zi) 2 , 



2 _ _ 2 



a, 2 < 



a 2 = {n + ry + (z + ZiY 



jpi is the toroidal component of the current in the center 
of the cell i, 

tt/2 



F(k, f ) = / —, ^^— 2 elliptic integral of the first 

q V (1 — fc 2 sin t) 

kind, 



tt/2 



LX\ K ,K, 2 ) - J (1 _ fe 2 sin 2 t) 3/2 

third kind. 

Taking into account, that J = j^rotH, and using (A. 5) 
we get boundary values of the poloidal magnetic field at 
the outer boundary. 

Approximation error for boundary values of the 
poloidal components of the magnetic field in the formu- 
lae ( |A.5|) is 0(h). However using these formulae requires 
0(Ny/ N) operations (where N is the total number of grid 
knots) and takes a lot of CPU time. During calculation the 
values of the poloidal components of magnetic field at the 
outer boundary change weakly and to decrease CPU time 



elliptic integral of the 



we used formula (A.5) every tenth time step. At all other 



time steps we extrapolate poloidal fields from the vicinity 
of the outer boundary to the outer boundary. 



U r0 = 0, u z0 = 0, po = 1, T = 0, 

H r0 = 2.507, H z0 = 1.401. 
The piston has a velocity: 

V r piston = 0.25, V z piston = 0.2735. 

The values of the parameters of the gas after the shock 
wave front (defined by the Hugoniot formulae) are the 
following: 

u r = 0.25, u z = -0.2735, T = 0.0117, p = 1.33, 

H r = 2.507, H z = 2.802. 

The discrepancy in the analytical and numerical values 
after the MHD shock was < 1% for the magnetic field and 
< 0.5% for the density. The same test has been made for 
the shock, moving along z-axis with similar discrepancy. 
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1 The numerical data, used in this test are taken from 



[gamarskii fc Popov (1992) 
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